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Abstract 

We  explore  the  praticability  of  optimal  shape  design  for  flows  modeled  by  the 
Euler  equations.  We  define  a  functional  whose  minimum  represents  the  optimality 
condition.  The  gradient  of  the  functional  with  respect  to  the  geometry  is  calculated 
with  the  Lagrange  multipliers,  which  are  determined  by  solving  a  costate  equation. 
The  optimization  problem  is  then  examined  by  comparing  the  performance  of  several 
gradient-based  optimization  algorithms.  In  this  formulation,  the  flow  field  can  be 
computed  to  an  arbitrary  order  of  accuracy.  Finally,  some  results  for  internal  flows 
with  embedded  shocks  are  presented,  including  a  case  for  which  the  solution  to  the 
inverse  problem  does  not  belong  to  the  design  space. 
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1  Introduction 


A  classical  problem  in  engineering  is  to  define  the  shape  of  a  vehicle  to  achieve  a  required 
performance.  In  fluid  dynamics,  techniques  have  been  developed  to  solve  the  following 
inverse  problem;  given  a  pressure  or  a  velocity  distribution  over  an  aerodynamic  body, 
determine  the  corresponding  geometry.  See,  for  example,  reference  [8].  A  broader  category 
of  problems  can  be  solved  by  means  of  optimization,  provided  that  one  is  ready  to  accept 
the  necessity  of  computing  the  flow  field  hundreds  of  times.  In  using  models  of  increased 
complexity  to  describe  the  flow  field,  the  development  of  new  algorithms  is  necessary  in 
many  cases  to  reduce  the  computational  load.  In  this  paper,  we  investigate  one  method  for 
achieving  this  reduction.  The  variational  technique  that  is  applied  in  this  work  has  been 
used  since  before  complex  flows  could  be  integrated  numerically.  See,  for  example,  refer¬ 
ence  [1].  Jameson  [4]  was  the  first  to  apply  this  technique  to  computational  fiuid  dynamics. 
With  this  approach,  a  functional  or  cost  function  is  determined  such  that  its  minimum 
represents  an  optimal  solution.  By  introducing  a  set  of  Lagrange  multipliers,  the  gradient 
of  the  functional  can  be  calculated  with  respect  to  the  geometry  by  computing  the  flow 
field  only  once  for  each  gradient  evaluation.  For  incompressible  irrotational  steady  flows,  a 
further  reduction  in  the  computational  effort  is  possible.  (See  [10].)  The  formulation  devel¬ 
oped  in  this  work  is  more  suited  to  a  correct  analysis.  In  fact,  this  formulation  eliminates 
the  need  to  consider  the  flow  field  variables  dependent  on  the  geometry.  In  the  present 
work,  we  extend  the  work  presented  in  reference  [5].  In  reference  [5],  an  exact  gradient 
of  the  functional  with  respect  to  the  design  variables  was  obtained  on  the  discrete  level, 
which  can  be  a  limitation  for  compressible  flows  because  in  presence  of  shocks  the  discrete 
functional  can  present  discontinuities.  In  considering  compressible  flows  with  embedded 
shocks,  we  derive  the  gradient  on  a  continuous  level.  Furthermore,  we  provide  a  method 
for  calculating  the  conditions  that  the  Lagrange  multipliers  must  satisfy  at  the  boundaries 
and  at  the  shock.  Finally,  we  point  out  that  our  formulation  can  be  used  with  complex 
flow  solvers  because  the  differentiability  of  the  solver  is  not  requested. 

2  Problem  Statement 

The  Euler  equations  are  given  by 


Ut  Fa;  -f  Gj,  —  0  (1) 
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u  =  X  component  of  velocity  vector 
V  —  y  component  of  velocity  vector 
e  =  specific  total  energy 
p  —  pressure 
a  =  speed  of  sound 
7  =  specific  heats  ratio 
7-1 


and  p  =  Kp{2e  -  -  v^).  Furthermore, 

F  =  -I^U  =  A(U)U  (2) 


and 


G  =  =  B(U)U 


(3) 


where  A  and  B  are  given  in  app.  I. 

We  assume  that  these  equations  are  defined  on  a  physical  space  $ .  In  this  space  is  included 
a  subdomain  whose  boundary  is  denoted  by  F.  On  the  boundary,  we  define  a  curvilinear 
coordinate  s  and  a  normal  n  =  {uxiriy)  that  points  outward.  (See  fig.  1.) 

The  optimization  problem  studied  here  is  defined  as  the  minimization  of  the  functional 
£  =  /p  (^(p,  p,  u,  u,  r)ds  over  all  admissible  shapes  of  the  subdomain  fl,  subject  to  the 
steady-state  Euler  equations  with  proper  boundary  conditions  on  F . 

Although  the  method  we  present  is  general,  we  focus  on  the  following  model  problem.  The 
subdomain  0  is  represented  by  a  nozzle.  (See  fig.  2).  At  the  inlet,  total  pressure,  total 
temperature,  and  the  ratio  <t  =  vju  are  fixed.  At  the  outlet,  if  the  flow  is  subsonic,  the 
static  pressure  is  fixed,  and  at  the  solid  walls  the  impermeability  condition  unx  -|-  vuy  = 
0  is  enforced.  The  upper  wall  is  kept  fixed.  The  lower  wall  0  is  represented  by  the 
parameterization 

!/(0)  =  E“i/<W  W 

i 


where  the  functions  fi{x)  are  some  shape  functions  and  o:  =  (ai, ..., a,, ...)  is  the  corre¬ 
sponding  set  of  shape  coefficients.  Given  a  desirable  lower  wall  pressure  distribution  p*(x) 
and  the  actual  pressure  distribution  on  the  lower  wall  p"'(a;),  the  optimization  problem 
consists  in  finding  a  set  of  shape  coefficients  Oj  such  that  the  functional 


(5) 


is  minimized. 
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3  Lagrange  Multipliers  and  Optimality 


The  problem  of  achieving  the  minimum  is  addressed  by  introducing  a  set  of  Lagrange 
multipliers.  Consider  the  augmented  functional 

£(U,  a,  K,  n)  =  E  I  ^A(AUa;  +  +  f  npY  -  nds  (6) 

JCl  J  Q 

where  V  =  (u,  u).  The  vector  A(x,  y)  =  *(Ai,  A2,  A3,  A4)  and  the  scalar  ij.{s)  are  the  contin¬ 
uous  equivalents  of  the  Lagrange  multipliers. 

We  calculate  the  variation  of  the  functional  C  with  respect  to  the  variation  of  the  functions 
U,  A,  and  p  and  the  parameters  ai,  respectively.  When  U(ar,?/)  is  increased  by  a  function 
eU(a;;  t/),  the  functional  C  increases  by  an  amount  e6Cu-  In  the  same  way,  A{x,y)  is 
increased  by  eA{x,y)\  p{s),  by  ep{s)]  and  each  Oj,  by  gd,-. 

If  we  follow  the  derivation  in  app.  II  and  take 


SjC  —  SC-jj  +  SEa.  “k  EjCu  T  d£,a 


then  we  obtain 


t  ^  iP"  -  P*)^dx  +  I  'A(An,  -b  Bny)Vds  + 
Ja  uVj  0 

j^(*A,A  +  ‘A,B)Udfi  +  J^pn^  TJds 
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dpV  _  /  0  1  0  0 
5U  “  \  0  0  1  0 


Furthermore, 


6Ca  =  [  ^A(AU^  +  BlJy)dn 

Jq 

SCy,  =  ppY  ■  nds 
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+  [  p  ^i^^  -nfids—  [  ppY  ’  cos^  0  ds 

Je  oy  Je  ax 

+  /  ///?V  •  n  ^  sin^cos^c?5  <5^ 

Je  ax 


where  6  is  the  angle  between  the  normal  n  and  the  y-axis  and  t  =  {—ny^rix)- 


At  the  minimum  of  the  functional,  we  have  =  0  for  all  possible  choices  of  the  functions 
U,  A,  and  fji  and  of  the  parameters  a.  This  condition  is  reached  when 

8Cu  =  S£a  =  6C^  =  =  0  (11) 

Note  that  because  of  the  necessary  conditions  (eqs.  (11))  the  unconstrained  minimum  of  the 
functional  £(U,a,  A,)w)  corresponds  to  the  constrained  minimum  of  the  functional  £{a). 
In  fact,  we  have 

6U=0  ^  AU,,  +  BUj,  =  0  and  =  0  ^  pV  •  n  =  0 

which  means  that  U  must  satisfy  the  Euler  equations  with  boundary  conditions.  In  addi¬ 
tion,  for  the  minimum  of  £  we  have  8£tj  =  0,  which  leads  to 

‘AA^  +  ^BAy  =  0  on  (12) 


^  (p"  -  p*)  cos  e  -b  *A(An^  -b  Bnj,)  -b  pn  =  0  on  0.  (13) 

At  the  inlet,  outlet,  and  upper  wall, 

'A(An^  -b  Bn,)U  =  0  (14) 

Given  U  and  the  set  of  costate  eqs.  (12)  and  (14),  we  can  determine  uniquely  A  in  and 
p  on  0.  (See  app.  III). 

Finally,  given  a  and  given  U  and  A  from  the  above  equations,  we  can  calculate  from  eq. 

(10) 


1  =  ^  (p-  -  p*)  fi  dx+  f  *A( AU,  +  BU,)  fi  cos  e  ds 

a  Ja  dy  0  Je 

+  /  P^^^^  ■  nfids  —  f  ppy  ■  t  cos^  6  ds 

Je  dy  Jq  dx 


f  dfi 

I  ppY  •  n  -jl-  sin  ^  cos  0  ds 
Je  dx 


In  cases  for  which  shock  occurs  in  the  flow  field,  we  split  the  domain  of  integration  by 
means  of  a  curve  T  that  coincides  with  the  shock  where  it  exists.  Then,  we  follow  the  same 
derivation  so  far  on  each  of  the  two  subdomains,  with  T  as  a  boundary.  (See  app.  III.) 
The  strategy  that  we  use  to  achieve  the  minimum  of  £  is  as  follows: 


1.  Start  with  a  set  a  of  shape  coefficients. 

2.  Enforce  8£a  =  0  and  =  0  by  finding  U  such  that  it  satisfies  the  steady-state 
Euler  equations  and  boundary  conditions. 

3.  Enforce  8£u  =  0  by  finding  A  such  that  it  satisfies  the  costate  equations  and 
boundary  conditions. 
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4.  Calculate  V„£.  If  =  0  then  we  have  deternained  the  minimum;  otherwise 
continue  to  steps  5  and  6. 

5.  Update  a  with  criteria  based  on  Vcj£. 

6.  Restart  from  step  2. 

4  Discrete  Problem 

We  introduce  a  discrete  grid  that  is  defined  as  (x;,  j/m)  =  {xo  +  ^Ax,t/(0)  +  mAj/),  where 
Ax  is  constant  and  Ay  is  a  constant  fraction  of  the  local  height  of  the  nozzle.  (See  fig.  3.) 
The  steady  solution  of  the  Euler  equations  is  obtained  with  a  time-dependent  technique, 
in  the  frame  of  an  explicit  finite- volume  code.  The  conservative  variables  U  are  computed 
at  the  cell  centers,  and  the  fluxes  F  and  G  are  evaluated  at  the  cell  interfaces  with  the 
approximate  Riemann  solver  in  reference  [7].  Second-order  accuracy  is  achieved  by  using 
an  essentially  nonoscillatory  scheme  [2].  The  flow-field  values  at  the  cell  interfaces,  used 
as  initial  conditions  for  the  Riemann  problem,  are  reconstructed  by  means  of  a  linear 
interpolation  and  a  minmod  limiter.  The  amplitude  of  the  integration  step  is  chosen  in 
accordance  with  the  Courant-Friedrichs-Lewy  (CFL)  condition. 

The  costate  equations  are  discretized  on  the  same  grid  presented  above.  Because  they  have 
no  conservative  form,  the  numerical  solution  is  obtained  with  a  finite- difference  scheme. 
We  introduce  a  set  of  curvilinear  coordinates  (/p(x,y)  and  ^(x,y).  The  costate  equations 
are  then  written 

UA,,  +  =  0  (16) 

where  A  =  and  B  =  AiJ^x  +  B^j,.  The  transformations  (f  and  are  defined  as 

{xi,  ym)^l  and  (x/, Vm)  ^  rn,  respectively. 

We  find  the  solution  of  eq.(16)  as  the  asymptotic  limit  of  a  time- dependent  technique. 
Consider  eq.  (16)  embedded  in  time  as 

±  Af  +  +  'BA^  =  0  (17) 

We  must  select  the  proper  sign  for  the  time  derivative.  The  inlet  and  outlet  boundary 
conditions  for  the  costate  equations  are  complementary  to  those  of  the  flow  field  equations, 
in  the  sense  that  if  the  number  of  boundary  conditions  for  the  flow  field  is  c,  then  the  number 
of  boundary  conditions  for  the  costate  equations  is  4— c.  Therefore,  the  above  equations  and 
boundary  conditions  are  well  posed  if  we  select  the  negative  sign  for  the  time  derivative. 
In  fact,  the  resulting  characteristic  pattern  is  mirror  symmetric  with  respect  to  that  of  the 
flow-field  equations. 

In  the  presence  of  a  shock  in  the  flow  field,  the  matrices  M  and  are  discontinuous.  In 
particular,  the  characteristic  pattern  at  the  shock  indicates  the  necessity  of  a  boundary 
condition  for  the  costate  equations.  For  further  discussion,  see  reference  [3]. 

The  costate  equations  are  linear  and  as  such  are  the  boundary  conditions.  We  exploit 
this  property  to  numerically  solve  these  equations.  Suppose  that  locally  we  separate  the 
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(18) 


variables  through  the  following  approximation: 


t)  =  A'(v?,  t)  4-  A"(V’,  t). 

This  separation  of  variables  means  that,  for  example,  in  a  Taylor  expansion  about  the  point 
{<-p,'4>)  we  disregard  all  terms  that  involve  a  cross  product  This  approximation  is  at 
least  first-order  accurate.  We  substitute  eq.  (18)  into  eq.  (17)  to  obtain 

-A',-A"*+UA;  +  ^eA;  =  o 

and  we  are  left  with  the  following  subproblems  in  one  dimension: 


—  A  t  +  MA^  —  0 

^  (19) 

1 

> 

+ 

> 

II 

O 

(20) 

Let  us  define  n<^  =  + +  V’D  and  +  Vy)- 

The  left  and  right  eigenvector  matrices  of  A  and  B  are  calculated  by  using  the  formulas  in 
App.  I  with  n  =  and  n  =  n^,  respectively.  After  eqs.  (19)  and  (20)  are  diagonalized, 
we  upwind  the  derivatives  of  the  characteristic  variables  according  to  the  signs  of  the 
corresponding  eigenvalues.  The  time  step  At  is  chosen  according  to  the  CFL  condition. 
This  method  can  be  regarded  as  a  two-dimensional  interpretation  of  the  method  presented 
in  reference  [6]. 

The  boundary  conditions  can  be  split  in  a  similar  way.  Consider,  for  example,  the  boundary 
condition  at  the  solid  wall.  Because  eq.  (19)  is  defined  along  the  wall,  the  characteristic 
variables  can  be  upwinded  according  to  the  corresponding  eigenvalues.  On  the  contrary, 
the  third  row  of  eq.  (20),  which  corresponds  to  the  characteristic  with  a  speed  of  -fa,  is 
replaced  by  the  boundary  condition  in  eq.  (AIII.5).  Note  that  the  contravariant  component 
of  the  speed  in  the  direction  ^  is  0;  therefore,  the  resulting  system  can  be  written  as 


A'W"  =  0 
Aw;  =  0 

rixA^X'o  +  TiyA^X'!  =  —{p^“  —  P*)  cos  0  —  n^jA^A^  —  UyA^X'^ 
A‘W;  =  aAt(J<Pl  +  xtl)A’"Wl 


(21) 


where  A(-)  is  the  forward  finite  increment  of  the  function  (•)  with  respect  to  the  super¬ 
scripted  variable  and 

/  AWi  \ 

Am 
V  AWT  / 


‘L„-^AA 


In  the  third  row  of  eq.  (21),  we  have  the  functions  of  A* A*,  which  are  computed  separately 
as  mentioned.  The  other  boundary  conditions  are  enforced  in  the  same  pattern  that  is 
presented  above. 


5  Optimization  Experiments 

The  optimization  problem  is  addressed  with  four  different  gradient-based  criteria. 

1.  Steepest  descent  (SDl).  The  shape  coefficients  are  updated  as  follows: 
cc,-  <—  ai  —  v{dCI doii),  where  u  is  a  given  parameter. 

2.  Steepest  descent  with  v  selected  as  shown  (SD2).  Because  we  know  the  gradient 

at  the  present  iteration,  we  can  use  a  tentative  v  and  can  compute  the 
gradient  Va£  .  By  calculating  VaC  •  VaC, ,  we  linearly  estimate  v  such  that 
eventually  VaC  •  Va£ '  =  0.  Each  step  of  the  optimization  requires  solution  of 
the  flow-field  and  costate  equations  twice. 

3.  The  BFGS  algorithm  (presented  in  reference  [9]).  This  algorithm  (BFGSl)  ac¬ 
counts  for  the  curvature  of  the  hypersurface  C.  The  shape  coefficients  are  up¬ 
dated  according  to  the  formula  a*-  e-  ai  —  udi,  where  d  =  (...,  dj-, ...)  is  the 
descent  direction  determined  by  d  =  H  and  H  is  an  estimate  of  the  inverse 
of  the  Hessian  of  the  functional. 

4.  The  BFGS  algorithm  (as  above)  with  a  linear  estimate  of  v  such  that  d  •  VaC"  = 
0.  Each  optimization  step  requires  solution  of  the  flow-field  and  costate  equations 
twice  (BFGS2). 


The  computations  are  performed  on  a  40  x  20  grid  unless  otherwise  specified.  Total  pressure 
and  total  temperature  at  the  inlet  are  taken  unitary  and  <t(0,  y)  =  0.  At  the  outlet,  the 
static  pressure  depends  on  the  test  case  considered.  For  the  lower  wall  ordinate  t/(0),  we 
have 


y(0)  = 


'  0 

^  TiUi 

i  0 


if  —0.5  <  a;  <  0 
if  0  <  a;  <  1 
if  1  <  a;  <  1.5 


The  optimization  consists  in  finding  the  four  shape  coefficients  ai  such  that  the  modulus  of 
the  gradient  is  0.  The  flow-field  and  costate  equations  are  iterated  until  the  residuals 
are  less  than  10“®. 

Consider  a  test  case  in  which  a  pressure  distribution  in  a  subsonic  nozzle  is  found  for  which 
the  outlet  pressure  is  0.9  of  the  inlet  total  pressure.  We  take  a  =  (2, 0,0,0)  and  define  the 
corresponding  configuration  as  the  optimal  configuration.  Then,  we  compute  the  flow  field 
and  determine  the  pressure  distribution  on  the  lower  wall.  This  pressure  distribution  p*  is 
the  one  we  want  to  recover  with  the  optimization  algorithm.  In  fig.  6,  the  target  flow  field 
and  the  starting  configuration,  obtained  with  a  =  (—2,0, 0,0),  are  shown  along  with  the 
convergence  histories  for  DSl  and  BFGS2. 

For  the  supersonic  case,  we  take  a  constant  section  channel  as  a  starting  configuration  and 
a  =  (2, 0,0, 0)  as  the  target.  In  fig.  7,  we  present  the  results  obtained  when  the  outlet 
pressure  is  lowered  to  0.5  of  the  inlet  total  pressure.  A  relevant  shock  is  generated  in  the 
flow  field.  In  the  first  optimization  iterations,  we  updated  the  shape  coefficients,  as  was 
done  for  DSl.  This  step  is  necessary  because  this  far  from  the  minimum  the  functional 


7 


C  might  be  not  convex;  therefore,  the  estimate  of  v  that  was  used  would  not  be  correct. 
Figure  8  shows  the  sequence  of  lower  wall  configurations  obtained  with  BFGS2.  We  also 
tested  the  effect  of  grid  coarseness  by  reducing  the  number  of  grid  points  to  20  x  10.  (See 
fig-  9.) 

The  second  test  case  is  designed  to  check  the  capability  of  the  algorithm  m  detecting 
minima  in  cases  for  which  the  desired  pressure  distribution  is  out  of  the  design  space  (i.e., 
the  functional  is  not  0  at  the  minimum).  The  pressure  distribution  f  is  obtained  with  an 
outlet  boundary  condition  that  differs  from  the  distribution  that  is  actually  used  in  the 
optimization  routine.  The  results  are  given  in  fig.  10. 

The  DSl  updating  strategy  had  the  least  attractive  rate  of  reduction  of  the  functional.  Our 
experience  showed,  nevertheless,  that  it  was  the  most  reliable  in  cases  of  complicated  surface 
topologies  that  can  occur  in  flow  fields  with  embedded  shocks.  The  BFGS2  becomes  the 
most  efficient  of  the  algorithms  tested  when  it  is  coupled  with  DSl.  With  this  algorithm, 
the  functional  was  reduced  by  orders  of  magnitude. 

Each  optimization  procedure  with  a  40  x  20  grid  required  6  h  of  central  processing  unit 
(cpu)  time  on  a  DEC  3000/500.  The  20  x  10  case  required  1  h  of  cpu  time. 


6  Concluding  Remarks 

We  have  derived  an  expression  of  the  gradient  of  the  cost  function  with  respect  to  the  shape 
coefficients.  The  boundary  conditions  for  the  costate  equations  have  been  presented;  we 
have  shown  their  relevance  to  the  well  posedness  of  the  problem.  In  the  case  of  shocks, 
we  provided  the  proper  conditions  for  the  costate  equations  at  the  discontinuity.  On  the 
discrete  level,  we  proposed  a  method  of  integrating  the  costate  equations  in  accordance  with 
a  revisited  scheme.  Additional  work  is  needed  to  test  the  algorithm  with  more  realistic  test 
cases  and  to  apply  the  One-Shot  method  (reference  [10])  to  hyperbolic  problems. 
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Appendix  I 

The  Jacobian  matrices  for  the  Euler  equation  in  conservative  variables  are 


B 


0 

kV^  - 
—uv 

(—76  +  2kV‘^)u  7e 

0 

—uv 
kV^  - 


1  0  0 

(3  —  7)u  —2kv  2k 

V  u  0 

-  kV^  —  2ku^  —2kuv  'ju  J 

0  10 

V  u  0 

—2ku  (3  —  7)u  2k 


(AM) 


(AI.2) 


[  (—76  +  2kV^)v  —2kuv  7e  —  kV^^  —  2kv^ 

The  Jacobian  in  the  direction  n  is  C  =  An^+Buy.  The  left  (L„)  and  right  (L'^)  eigenvector 
matrices  of  C  are 


1  - 

kV'‘la? 

2kul  a? 

2kv/a^  —2kla?  1 

v.Ip 

Uylp 

0 

1 

(-K.  +  kVya)/p 

(n^  -  2kula)/p 

(ny  —  2kvla)fp  2klpa 

(V„  +  kVya)lp 

-{n^  +  2kuja)fp  - 

-[uy  +  2kv !  a) !  p  2klpa 

-  1 

0 

p/2a 

pl2a 

u 

pUy 

p{u  +  anx)l2a 

p{u  —  anx)l2a 

V 

prix 

p{v  +  any)f2a 

piy  —  any)l2a 

.  vy2 

-pVt  p{V^  +  aVA:  +  2aVn)/ia 

p{V^  +  a^k  -  2aK)/4c  . 

where  I4  =  V  •  n  and  y,  =  V  •  t.  The  diagonal  matrix  Dn  =  L„CnL„^  is 

'Vn  0  0  0  ■ 

OK  0  0 

0  0  K  +  «  0 

0  0  0  K  - «  . 


Appendix  II 

To  calculate  SCu,  consider  the  increment  U  ^  U  +  eU 

F  < —  F  +  sAU  +  h.o.t.  and  G  ■« —  G  +  sBU  +  h.o.t. 


(AI.3) 


(AI.4) 


(AI.5) 


We  obtain 


SCu  =  / 

Ja 

L 


dV 


e 
dpY 

e  oV 


{p  -  p*)Udx  +  /  'A  [(AU)^  +  BU)^J  da 

J  Q 

\Jds  +  h.o.t. 


(A11.1) 
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If  we  apply  Gauss’  theorem  to  the  second  integral  of  the  above  equation  then  we  find  eq. 
(7).  Equations  (8)  and  (9)  are  easily  obtained. 

To  calculate  SCa,  we  first  consider  the  variation  of  the  functions  defined  on  0: 


Then,  we  consider  the  variation  of  the  geometry;  other  higher  order  effects  are  disregarded. 
When  the  geometry  is  perturbed,  the  domain  of  integration  fi,  the  normal  n  and  the 
element  of  integration  ds  are  perturbed.  The  domain  is  increased  (fig.  4)  by  a  quan¬ 
tity  eoLifi  cos  Ods.  The  normal  n  is  perturbed  by  a  quantity  -s&idfildx  cos^^t;  ds,  by 
eaidfifdx  cos  6  smO  ds.  (See  fig.  5.) 


Appendix  III 

Consider  eq.  (14).  This  elation  defines  the  boundary  conditions  for  A  after  we  impose 
the  proper  constraints  on  U.  At  the  inlet,  only  one  compon^t  of  the  variation  of  the  flux 
in  the  direction  normal  to  the  boundary  Fn  =  {An^  +  is  independent  of  the  others 

because  total  pressure,  total  temperature,  and  a  are  fixed.  If  we  express  all  components  of 
F„  in  terms  of  we  obtain 

[C  +  -  AI^)] 

K-7/(1 

m 

where  ^  -|-  any,  rj  =  {uy  -  auj;),  H  is  the  specific  total  enthalpy,  and  M  is  the  local 

Mach  number.  For  an  arbitrary  choice  of  from  eq.  (14)  we  have 

Ai^  +  X2U  +  w/i'i-  -  M^)]  +  Xzu  [aC  -  v/il  -  M^)]  +  X^Hi  =  0  (AIII.l) 

At  the  outlet,  if  the  flow  i^^subsonic,  then  only  the  static  pressure  is  fixed;  therefore,  three 
components  of  the  vector  U  are  arbitrary.  If  we  take  ^  as  the  dependent  variable,  we  have 

punj,  -1-  '^Uy 

pu(]4  -1-  UUx)  +  UUy'^  -  puVn 

^{Vn  -h  VUy)  -1-  -  pvVn 

p{'ypf2kp  +  +  V^)  -f  ^{VnU  +  HUx)  +  ^{VnV  Huy) 

For  an  arbitrary  choice  of  p,  pu,  and  pv,  from  eq.  (14)  we  have 

XiUx  H"  A2(V^  d"  un-x)  +  X^vUx  A4(I4i'w  +  Hnx)  —  0 
XlUy  X2Uny  +  A3(y„  +  vny)  -f  A4(K^^  HUy)  =  0  (AIII.2) 

A2UK  +  A3UK  +  XiVn{ipl2kp  -\-u^  -\-v^)  =  0 
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For  a  supersonic  outlet,  no  conditions  exist  on  U;  therefore,  we,  have 

A  =  0. 


(AIII.3) 


If  a  shock  is  embedded  in  the  flow  field,  then  the  shock  is  considered  as  a  boundary  for 
the  costate  equations.  The  consequent  boundary  conditions  are  applied  on  each  side  of  eq. 
(AIIL3)  before  the  shock  and  eq.  on  each  side  of  (AIII.l)  after  the  shock. 

For  the  upper  wall,  we  have 

r  ®  ' 

Fn  =  ..  P 


such  that  for  an  arbitrary  choice  of  p  eq.  (14)  is  satisfied  if 

^2'^x  T  —  0 


(AIII.4) 


At  the  lower  wall,  eq.  (13)  applies.  Because  the  rank  of  Aux  +  Buj,  at  the  wall  is  2,  the 
system  has  only  two  linearly  independent  rows.  We  obtain 


+  (p®  -  p*)  cos$  =0 


(AIII.5) 


which  is  the  boundary  condition  for  A,  and 

fi  =  —  Ai  +  uX2  +  vXs  +  {■ye  —  ^y^)A4j  (AIII.6) 

which  is  the  relation  between  p  and  A  on  the  boundary. 


12 


Figure  1:  Physical  space. 


Figure  2:  Model  problem. 


(a)  Target  Mach-number  field. 


(b)  Starting  configuration. 


(c)  Functional  and  modulus  of  gradient 
versus  number  of  iterations  for  SDl. 


(d)  Functional  and  modulus  of  gradient 
versus  number  of  iterations  for  BFGS2. 


,5th  Iter 
JOth  Iter 
.15th  Iter 
<20th  Iter 
>25th  Iter 
^30th  Iter 
535th  Iter 
z40th  Iter 
;;45th  Iter 


Figure  8:  Wall  shapes  sequence  with  BFGS2.  Starting  configuration;  constant  section. 
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